FAST SIMULATION OF TRUNCATED GAUSSIAN 
DISTRIBUTIONS 



NICOLAS CHOPIN, ENSAE-CREST 



Abstract. We consider the problem of simulating a Gaussian vector X, con- 
ditional on the fact that each component of X belongs to a finite interval 
[cti,bi], or a semi-finite interval [oj,+oo). In the one-dimensional case, we 
design a table-based algorithm that is computationally faster than alternative 
algorithms. In the two-dimensional case, we design an accept-reject algorithm. 
According to our calculations and our numerical studies, the acceptance rate 
of this algorithm is bounded from below by 0.5 for semi-finite truncation in- 
tervals, and by 0.47 for finite intervals. Extension to 3 or more dimensions is 
discussed. 



1. Introduction 



Let X = (Xi, . . . ,Xd) be a d— dimensional Gaussian vector with mean ii and 
covariance matrix E, and let [a,i,bi] be d intervals, where hi may be either a real 
number or +oo. The distribution of X, conditional on the event th at Xj € \a,j, bj , 
i = 1, ...,d, is usually called a truncated Gaussian distribution ([Johnson et al 
13) 



0, and that 



1994 Chap. 13). Without loss of generality, one may assume that \i 
E has unit diagonal elements. 

Numerous statistical algorithms rely on intensive simulation of truncated Gaus- 
sian distributions. In particular, several Bayesian models generate full conditional 
distri butions of this type, either directly or through a data augmentation represent- 
ation ( Tanner and Wong . 1987t ) . Thus, the corresponding Gibbs samplers (or more 
generally Markov chain Monte Carlo algorithms) draw repetitively from truncated 
Gaussian di stributions . Examples include linear regression model s with ordered 
param eters ( Chen and Deelvj. 1996) or appl ied to truncated data ( Gelfand et al. 



1992 ) , probit models ( Albert and Chibl. 19931) . m ultinomial probit models ( Albert and Chibl. 

1993 ; McCulloch and Rossi , 1994 ; NobileL 199S ), multivariate probit m odels ( Chib and Greenbergl 
1998 ), multiran ked probit models ( Linardakis and Dellaportasl. 12003]) . to bit mod- 
els (jChib . 1992 ). models used in sp ectroscopy ( Gulam Razul et al. . 2003 ). copula 



regression models (jPitt et all 120061 ), among others 



To understand how intensive such MCMC algorithms can be, consider the prob- 
lem of sampling the posterior distribution of a multinomial probit model with n 
observatio ns and p alternatives. A solut ion is to perform T iterations of the Gibbs 
sampler of McCulloch and Rossil (1994), but this requires the generation of Trip 



univariate truncated Gaussian variates, a number that may exceed 10 or even 
10 15 in difficult scenarios. Hence any improvement with respect to the computa- 
tional cost of simulating univariate truncated Gaussian distributions may lead to 
important savings. Another important aspect of such algorithms is that they simu- 
late only one random variable from a given truncated Gaussian distribution, that is, 
the parameters and the truncation intervals [aj,6j] change every time a truncated 
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Gaussian variate is generated. Thus, we are interested in developing specialised 
algorithms which are guaranteed to generate quickly one random variable from the 
desired distribution, for all possible inputs (i.e., parameters and truncation inter- 
vals). As a corollary, these algorithms cannot afford a long set-up (initialisation) 
time, where some exploration of the target density is performed in order to improve 
performance; this type of initialisation is meaningful only when one needs to simu- 
late many variables from one fixed distribution, and is not discussed in this paper. 
The algorithm we propose does require a table set-up, but which is independent of 
the input parameters. 

The first part of this paper presents a table-based simulation algorithm for uni- 
variate Gaussian distributions truncated to either a finite interval [a, b] or a semi- 
finite interval [a, +00). In the latter case, and given the truncation point a, our 
algorithm is up to three times faster than alternative algorithms in our simulations; 
se e below for references . Our algori thm is inspired from the Ziggurat algorithm 
of Marsaglia and Tsane ( 1984 . 200dh . whi ch is usually c onsidered as the fastest 
Gaussian sampler, and is also very close to Ahrens ( 1995h 's algorithm. 

Another possible strategy for accelerating a Gibbs sampler is to 'block', i.e., to 
update jointly, two or more components of the posterior density; this often strongly 
improves the mixing properties of the algorithm. In some of the aforementioned 
models, blocking requires simulating multivariate truncated Gaussian variates. We 
develop an accept-reject algorithm for simulating from bivariate truncated Gaussian 
distributions. In all but one particular case for finite intervals, we manage to prove 
formally that the acceptance rate is bounded from below by 0.22. Our numerical 
studies seem to indicate that this bound is not optimal, and that the acceptance rate 
is bounded from below by 1/2 when the truncation intervals are semi-finite, and by 
0.477 when they are finite. (This remains true even when the correlation coefficient 
get close to 1 or — 1, that is, in situations where MCMC blocking is particularly 
efficient.) In the former case, we explain how to generalise this algorithm in some 
situations to truncated Gaussian distributions of dimension d, with the outcome 
that the acceptance rate is bounded from below by 1/2 . Interestingly, some of 
the constants that must be pre-computed for our univariate algorithm can be re- 
used so as to bypass part of the computations performed by our multi-dimensional 
algorithms. 

We note that independent variables from trun cated Gaussian distributions may 
also be obtained using the perfect samplers of iPhilippe and Robert! <|2003h and 



Fernandez et al.l (|2007f) . but, for small dimensions, these algorithms are much more 



expensive than our approach, since each sa mple requires running a Markov chain un- 

til som e criterion is fulfilled. (According to Hormann and Levdoldl ( 2006f ) . Philippe and Robert 
(|2003l ) may not sample from the correct distribution.) 

The paper is organised as follows. Section 2 presents our algorithm for simulating 
univariate truncated Gaussian variables. Section 3 presents a rejection algorithm 
for simulating bivariate Gaussian vectors, the components of which are truncated 
to semi-finite intervals [a*, +00). Section 4 does the same thing for finite truncation 
intervals. Section 5 explains how to generalise the algorithms of Section 3 to three 
or more dimensions. Section 6 concludes. 
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2. One-dimensional case 



First, we consider the problem of simulating a random variable X from a uni- 
variate Gaussian density truncated to [a, +00): 



(2.1) 



p(x) 



(p(x) 



I(x > a) 



$(-a) 

for some truncation point a, where if and $ denote respectively the unit Gaussian 
probability density and cumulative distribution functions; <p(x) — exp(— x 2 /2)/^/2tt. 
The extension to a finite truncation interval [a, b] is explained in H2.5I 



2.1. Review of current algorithms. A convenient way to generate X is to use 
the inverse transform method: 

(2.2) X = (<S>(-a)U) , 

where U ~ U[0, 1] is a uniform variate. Note that this expression is equivalent to 

(2.3) X = $- 1 ($(a) + {l-$(a)}[/), 

but the latter expression is less stable numerically for large values of a, because it is 
easier to approximate $~ 1 in the left tail than in the right tail. In our experiments, 
(|2.3|) generates "inf" values when a > 9.5, while (|2.2j) generates "inf" values only 



when a > 37.5. 

As noted by iGlassermanl (|2004 Chap. 2), the inverse transform method seldom 
produces the fastest algorithms, but it has appealing properties that may justify 
the increased cost in some settings, in particular when used in conjunction with 
variance reduction o r quas i Monte Carlo techniques; see the same reference and 
also e.g. Blair et al. ( 19761 ) for an overview of fast methods for evaluating $ and 
We now focus on specialised algorithms. 

First, we recall br iefly the rejection principle (e.g. Devrovel 1 19861 Chap. 2 or 
Hormann et al. J2004 Chap. 2). Assume we know of a proposal density q such that 



p(x) < Mq(x) 

for some M > 1, and all x in the support of q. Then a sample from p can be 
obtained as follows: simulate X ~ q, and accept the realisation x with probability 
p(x)/Mq(x): otherwise repeat. The expected acceptance probability, a.k.a. the 
acceptance rate, equals 1/M. It is important to choose q so that a) M is small and 
b) simulating from q is cheap . 

For a > 0, Devrovel (1986, p. 382) proposes a rejection algorithm based on the 
proposal exponential density q(x) = A exp {— A(x — a)}, for x > a, with A = a. 
The acceptance rate of this algorithm is aexp(a 2 /2)$(a), which goes to zero as 
a — > 0, so it can be used only for a > clq, with say oo = 1. For a < do, one 
may use i nstead the follo wing trivial rejectio n algorithm : repe at X ~ N(0, 1) until 



X > a. Devrove ( 19861 p. 382) mentions Marsaglia ( 1964 ) 's algori t hm, w hich 



has the same acceptance rate, but is a bit more expensive. Gewekd (1991) and 
iRoberd (|l995h independently derive a rejection algorithm for a > 0, based again on 
q(x) = Aexp{— X(x — a)}, but with A = (a + \J a 2 + 4)/2, which is shown to give 
the optimal acceptance rate. For a < 0, these authors use the sam e trivial sampler 
as ab ove. In pr inciple, th e se alg orith ms may be refined using A RS (jGilks and Wild, 
1992 ). see also Hormann ( 1995f l and lEvans and Swartz ( 19981 ). i.e., rejected points 



are used to improve the proposal density (which is then piecewise exponential) . As 
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FIGURE 2.1. Plot of N(0, 1) density and the 2N vertical rect- 
angles, for N = 20. 



explained in the introduction however, we are interested in situations when only 
one random variate must be generated (for a given value of a); hence, since the 
acceptance rate of Devroye's and Geweke and Robert's algorithms are high enough, 
we do not discuss this further. 

The algorithm we propose in this paper is faster than these specialised algorithms 
for two reasons: (a) its acceptance rate is higher, and, in fact, is almost one, for 
most values of a; and (b) with high probability, the only floating point operations 
that the algorithm performs are 2 additions and 3 multiplications, whereas the 
aforementioned algorithms computes a few logarithms and square roots. 

2.2. Principle of proposed algorithm. For the sake of clarity, we consider first 
the simulation of a non-truncated N(0, 1) density, and consider the extension to a 
truncated density in next section. We do not claim, however, that this algorithm is 
either interesting or novel in the non-truncated case, see below for references. The 
principle of the algorithm is summarised by Figure [2~T1 The proposal distribution 
consists of 2N + 2 regions: 2N vertical rectangles of equal area, and two Gaussian 
tails of the same area. For rectangle i, i — —N, . . . , N — 1, let [xi, Xi+i] denote its 
left and right x-ordinates, yi its height, i.e., yi — <p(xi) V cp(xi+i), y. the height 
of the smaller of its two immediate neighbours, i.e., y. = tp(xi) A ip(x i+ i), and let 
di = Xi + \ — Xi, Si = diyi/y^ (Symbols A and V means 'min' and max' throughout 
the paper.) All these numbers are computed beforehand and defined as constants 
in the program. Note that the region labelled — N — 1 (resp. N) is the left tail 
(resp. right tail) truncated at x = X-n (resp. at x = xn). 

To sample X ~ N(0, 1), one may proceed as follows: choose randomly region 
i, sample the point (X, Y) uniformly within the chosen region, and accept X if 
Y < (f(X); otherwise repeat. However, if the chosen region is a rectangle, most of 
the computation can be bypassed: one may first simulate Y, i.e., draw U ~ U[0, 1] 
and set Y — yJJ , without simulating X, and check that the realisation y of Y is 
such that y < y:, recall that y . = tp(xi) A ip(xi+i). If this condition is fulfilled, then 
the realised pair (x, y) must be accepted whatever the value of x. Furthermore, one 
can recycle the realisation « of [/, and therefore avoid drawing a second uniform 
variate, by simply setting x — Xi + 5iU. 

In short, with high probability, the algorithm only performs the following basic 
operations: 
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draw a random integer i uniformly in range —N— 1, . . . , N 
if i < N and i > -N - 1 then 

sample u ~ U[0, 1] 

y <- iji * u 

if y < 2/^ then 

return x.j + S t * u 

end if 
end if 



A complete outline of the algorithm is given in Appendix A. When the condition 
y < y . is not fulfilled, one must sample X and check that (X, Y) is indeed under 
the curve of the N(0, 1) density. Similarly, when the chosen region is either the left 
or right tail, one may use Devroye's algorithm in order to simulate X. 

Note that this algorithm has a slightly higher numerical precision than the rejec- 
tion algorithms mentioned in the previous section: when calculating x — Xi +di *w, 
the absolute error equals the precision of the random generator that produces u, 
say 2 -32 « 2.3 x 10~ 10 for a 32 bit generator, times the small number Si, which is 
typically of order 10 -3 . 

Many algorithms proposed i n the literature already use h i stogram s to co n struct a 
good p roposal density; see e.glM arsagli a and Tsand (jl984h . lAhrensI (jl993l ). IZamanI 



(jl996l ) or the survey in iHormann et alJ (|2004l . Chap. 5). In particular, the abov e 
algorithm is similar to the Ziggurat algorithm (jMarsaglia and Tsan 2 Ii984 I2000L 



which is the default Gaussian sampler in much mathe matical s oftwa re, e.g. Mat- 
lab or the GNU Scientific Library, and most similar to I Ahrend (|1995l )'s algorithm. 
Both algorithms already use the idea of using rectangles of equal areas, but the Zig- 
gurat algorithm is based on horizontal rectangles, while Ahrensl ( 19951 ) 's algorithm 



is based on vertical ones, as above. This seemingly innocuous variation greatly 
facilitates the extension to truncated densities, as explained in next section. 

2.3. Extension to truncated Gaussians. For a fixed truncation point a, let l a 
denote the index of the region that contains a. To adapt the above algorithm to 
the truncated density (|2.ip . one may choose an integer i a such that i a < l a , sample 
randomly one region among i a ,i a + 1, • • • , N, and proceed as explained above. In 
addition, for the regions i a ,...,l a , one must reject the random point [X, Y) if 
X < a, as described in Appendix A. 

The difficulty is to define i a in such a way that (a) the computation of i a is 
quick, and (b) i a is as close as possible to l ai so that the overall acceptance rate is 
as high as possible. We propose the following method. We choose a small width 
h > 0, and store beforehand in an integer array the following quantities: 

jk = max{i : xi < kh} , for all k such that kh 6 [a m in>«max] 
for some interval [a m i n ,a max ]. As said before, all these constants are computed 
separately, and hard-coded in the program. Then, provided a £ [a m i n ,a max ], i a is 
computed as 

ia = jy a /h\ 

where |_-J stands for the floor function. Provided h < min^Xi+i — x.i), each interval 
[kh, (k + l)h) contains at most one Xi, so that either i a = l a or i a — l a — 1. This 
means that, when choosing randomly between regions i a , i a + 1, . . . , N, one must 
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FIGURE 2.2. Execution time (seconds) vs truncation point, for 10 s 
simulations of Devroye's algorithm (dotted line), Geweke-Robert 
sampler (dashed line), the inverse tran sform alg o rithm (based on 
the inverse transform approximation of IWichura (1988), as imple- 
mented in the GSL library, dash-dotted line), and our algorithm 
(solid lines), with, from left to right, N s = 1000, N s = 2000, 
N. = 4000. 



treat separately the two leftmost regions i a and i a + 1, and perform the additional 
check mentioned above, i.e., x > a, but the other regions i a + 2, ...,N can be 
treated exactly as explained in the previous section. 

In our simulations, we set a m ; n = —2, a max = xn-20, and h = x\ — xq, that is, 
the smallest of the interval ranges (xj+i — Xi). A complete outline of the algorithm 
is given in Appendix A. 

2.4. Results. We implemented our algorithm, the inverse transform algorithm, 
Devroye's algorithm (using cut-off value ao = 0.65) and Geweke-Robert 's algorithm 
in C, using the GNU Scientific library (GSL). Figure l2~2l plots the execution time of 
10 8 runs on a 2.8 Ghz desktop computer, for each algorithm and for different values 
of the truncation point a. Our algorithm appears to be up to two times faster than 
Geweke-Robert 's algorithm, and up to three times faster than Devroye's algorithm 
and the inverse transform method. The three solid lines correspond to different 
sizes N s , N s — 1000, 2000, 4000, from left to right, of the five arrays containing 
the constants Xi, yi, y., d%, Sf, note that only those values such that Xi,Xi+i G 



[a m i n ,a max ] need to be stored, hence N s < 2N. Figure I2T2I shows that increasing 
N s only improves the execution time for a tiny interval of a values, so there may be 
little point in increasing N s further than, say, 4000. For N s = 4000, the acceptance 
rate is typically above 0.99 or even 0.999 for most values of a € [a m i n , dmax]- The 
increase of the computational cost for a > 2 is due to the increasing probability 
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of sampling X from the right tail using Devroye's algorithm. Outside of the range 
of the represented interval, the rejection algorithms have a similar computational 
cost, as they perform the same operations. For N s = 4000 and a double precision 
implementation, the total memory cost of the algorithm is 162 kB, a small fraction 
of the memory cache of most modern CPU's. (A CPU memory cache is a small, 
fast memory where a CPU stores data used repetitively.) 

2.5. Truncation to a finite interval [a, b]. The extension to finite truncation 
intervals [a, b] is straightforward. First, one determines a region index i a (resp. i' b ) 
such that either region i a or region i a + 1 (resp. region i' b or region i' b — 1) contains 
a (resp. b). using the table look-up method described in 12.31 Then, one proceeds 
as above, choosing randomly a region in the range i a , . . . , ii, and so on. 

However, a difficulty arises if (b — a) is small. Suppose for instance that a and 
b fall in the same region, and that b — a is small with respect to the width of 
the region. Then, if one samples uniformly point [X, Y) within that region, the 
probability that a < X < b may be arbitrarily small. 

We propose the following work-aroundi when — z a ^ ^mini say &min — 5, use 
instead a rejection algorithm based on iDevrovel (|l986lV s exponential proposal, but 
truncated to [a,b], i.e., 

(2.4) q (x) = _^£hM J(a < x < b) . 

exp(— Xa) — exp(— Xbj 

with X = a (resp. X = b) when b > (resp. when b < 0). The advantage of this 
approach is that it gives an acceptance rate close to 1 whatever the values of a, 
and 6, subject to i' b — i a < k m i n ; i.e., whether a and b are both in the same tail, 
or both close to 0. In the latter case, q should be close numerically to a uniform 
distribution. 

3. Bl-DIMENSIONAL CASE: SEMI-FINITE INTERVALS 

We now consider the simulation of X = (X\, X2) ~ ^(a 1 ; subject to X\ > a\ 
and X2 > ci2- Without loss of generality, we set fi — (0, 0)', 

1 P 
P 1 

and assume that a\ > ai\ if necessary, swap components to impose the last condi- 
tion. The joint density of the considered truncated density is, up to a constant: 

(3.1) v{x\,x%) oc exp | _ 2~2 ( x i + x 2~ 2p x i x 2) \ 1 ( x i > ai; x 2 > a 2 ) , 

where the short-hand v 2 = 1 — p 2 will be used throughout the rest of the paper. 
The conditional distribution of X 2 \X\ = x\ is a univariate Gaussian N(pxi,v 2 ) 
truncated to X 2 > a 2l which we denote from now on TNi a2 ^(pxi, v 2 ). A common 
misconception is that the marginal density of X\ is also a truncated Gaussian 
density, although standard calculus leads to: 

(3.2) p( Xl ) oc v ( Xl )$ ( PX1 ~ a2 ) I( Xl > a t ). 



E = 



In order to simulate from Q3.ip , a natural strategy is to derive a rejection al- 
gorithm for the marginal (|3.2p . and to simulate X 2 conditional on X\, using the 
algorithm we developed in Section 2. This is basically the approach adopted here, 
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although we shall see that, in some cases, it is preferable to derive a rejection 
sampler for the joint di stribution 1 3.111. W e mention briefl y that universal bivariate 
samplers exist, see e.g. IHormannl (|2000l) or iLevdoldl [20001. but as in the univariate 



case our objective is to design a specialised algorithm that runs faster (i.e., does 
not require a set-up time), for situations where only one random vector must be 
generated. 

To derive a proposal distribution for (|3.2| . we substitute the <!>(•) factor with a 
simpler expression derived from the two following straightforward inequalities: 

(3.3) i < <f>(x) < 1 for x > 0, 

(3.4) $(x) < c(x )ip(x) for x < x < 0, 

where c(xq) — A (— l/xo), for xq < 0, c(0) = We now distinguish 

between cases where the argument of $(•) in f|3 . 2(1 is positive, negative, or both, 
over the range of possible values for x%. We consider the following cases, and treat 
them separately: 

• case S + : either p > and pa\ — a 2 > 0, or p < and a\ < $ _1 (l/3) ss 
-0.4307 

• case S~: p<0, pa x - a 2 < 0, and ai > ^(l/S) « -0.4307. 

• case M + : p > and pa\ — a 2 < 0. 

• case M~: p < 0, pa x - a 2 > 0, and a x > $- 1 (-l/3) -0.4307. 
where 'S' stands for 'Simple', and 'M' for 'Mixture', as we elaborate below. 

We now prove that, in each case, it is possible to derive a rejection algorithm 
with an acceptance rate bounded from below for all values of p, a\ and a 2 . 

3.1. Case S + . Assuming first p > and pa\ — a 2 > 0, then, according to (|3.3p . 

(3.5) *(^l^ &[l /2,l] 

for all x > 01, which suggests the following proposal distribution: 

qs+{x\) oc ip{x\)I{x\ > Oi), 

i.e., a TiV[ ai _ +oo )(0, 1) distribution, in order to sample from the marginal p{x\). For 
a given x\ simulated from 95+, the acceptance probability equals (|3.5[) . hence the 
acceptance rate of such a rejection algorithm equals 

(3 - 6) y Q1 J^) dxi ' 

and is larger than 1/2 by construction. 

However, it is more efficient to simulate jointly (Xi,X 2 ) as follows: sample 
X\ ~ riY"r ai)OO )(0, 1), = jci ~ N(pxi,u 2 ), and accept if X 2 > 02; otherwise 

repeat. It is easy to check that the latter rejection algorithm has exactly the same 
acceptance rate, i.e., (|3.6p , as the former, but it is faster, because it does not perform 
any evaluation of $, and because X 2 is obtained 'for free', i.e., a second step is not 
required to generate X 2 . 

When p < 0, the argument of $ in (|3.6[) is not positive for all x% > ax, but 
the integral is still larger than 1/2 provided a\ < $ _1 (l/3). To establish this 
property, one may remark that (|3.6j) is the probability that X 2 > a 2 , conditional 
on X\ > 01, provided (Xi,X 2 ) ~ N 2 (p, E). This probability decreases with respect 
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to <22, 02 < 01, and, for 02 =01, this probability decreases with respect to —p and 
ai. Finally, for p = — 1 and ai = 02 = $ _1 (l/3), this probability equals 1/2. Thus, 
we use Algorithm S + also when p < and ai < $ _1 (l/3). 

3.2. Case S~ . If p < and pa\ — 0,2 < 0, then inequality (|3.4p holds for all values 
of the argument x — {px\ — 0,2) /v, and for xq = {pa\ — 0,2) /v- This suggests the 
following algorithm: sample X\ from proposal density 

q s -{xi) oc <p(xi)<p[— —)I(xi>ai) 



V v 

oc ip {x\] pa 2 ,v 2 ) I{xi> ai) 

that is, the density of truncated Gaussian distribution TNt aij oa)(P a 2, v 2 ), and ac- 
cept with probability: 

(3.7) ^(_P^~a2\ jc ( P a,-a2 



v 

where ip( x ) = x)/ip(x). The acceptance rate is then 



(3-8) ^TAT [aii00) ( / 9a 2 ,^ 2 ) 



We show formally in Appendix Bl that this acceptance rate admits a lower bound 
which is larger than 0.416; our numerical studies indicate that the optimal lower 
bound may be 1/2. Intuitively, the idea behind inequality (|3.4p is that $(— x) ps 
ip(x)/x for large values of x, hence the true marginal density p(x\) behaves like a 
Gaussian density density times 1/(02 — pX\), but the latter factor varies slowly with 
respect to a Gaussian density, so it can be discarded in the proposal. 

3.3. Case ill". If p < and pa± — 0,2 > 0, the quantity {px\ — a%)jv takes positive 
and negative values for X\ > a±. To combine both inequalities, one may use a 
mixture proposal: 
(3.9) 

q M - (xt) oc <p(xi)I (px 1 - a 2 > 0) + ^J^(p(xi)<p a2 ^j j ( pXl _ fl2 < 0) 

subject to x\ > a\. To sample from the mixture proposal, choose component 1 
(corresponding to the first term above), with probability lo\/ [uj\ + u> 2 )i with 

uji = <t>(a 2 /p) - $(ai) 

v ( a 2 , 1 (0,2V 

and choose component 2 otherwise. If component 1 is chosen, one can use the same 
shortcut as in Algorithm S + , that is, draw X\ ~ T~Ny ai M2 / p ] (0, 1) and X2IX1 = 
xi ~ N(px\, v 2 ), and accept the simulated pair (xi, X2) if 2:2 > °2- If component 2 is 
chosen, the proposed value for X\ is drawn from a TN\a 2 / Pt00 \{pa 2 , v ) distribution, 
and is accepted with probability 

£ , / pxi - a 2 

7T \ 

When a draw xi is accepted, it is completed with X2 drawn from X 2 \Xi — x± 
TN [a2 . +oo) (pxi,v 2 ). 



FAST SIMULATION OF TRUNCATED GAUSSIAN DISTRIBUTIONS 



10 



The acceptance rate of this algorithm is a weighted average (with weights given 
by uj\ and u) 2 ) of the acceptance rate of Component 1, which is larger than 1/2 by 
construction, and the acceptance rate of algorithm S~ for a 2 — pax, which is also 
bounded from below, as explained in the previous subsection. 



3.4. Case M + . If p > and pa\ — a 2 < 0, then again (px\ — a 2 )/v take both 
negative and positive values for x\ > ax, which suggests that we use a mixture 
proposal similar to (|3.9II . Unfortunately, the acceptance rate may be arbitrarily 
small in that case. Exact calculations are omitted for the sake of space, but it can 
be shown that the mode of p{x±) can be arbitrary far from x\ = a 2 j ' p, the point 
where the two components intersect, which gives an arbitrary small acceptance rate. 
Instead, we substitute p.4[) with a slightly different inequality: 

(3.10) §{x) < d{x )(p(x)e Xx for x < x < 

where A = 0.68, d(x ) = (^/nj2) V x(-x ), and = e Xx $(-x)/tp(x). This 

inequality stems from straightforward calculus. Other values of A are also valid, 
but in our numerical experiments, A = 0.68 seemed to be close to optimal, in terms 
of minimum acceptance rate. 

This inequality leads to the following proposal mixture density: 

Qm+(xi) « <p{n\)I (p%i - a 2 > 0) 

(px 1 -a 2 \ f X(pxi - a 2 )\ , ( pa\ - a 2 \ T . nN 
+ip(x\)ip I 1 exp I I d I I I (pxi - a 2 < 0) 

subject to xi > ax- The second term is proportional to a N(9,v 2 ) density, with 
9 = p(a 2 + \i>). To sample from this mixture, choose component 1, with probability 
ti/{t\ + t 2 ), choose component 2 otherwise, where 

n = $(-a 2 /p) 



T2 = "U(2t£zl\-9f*L 



'2tt { \ v ) \ v 

( 9 2 — a\ — 2\va 2 1 / pai — a 2 
CXP \ 2^ \ d {—^- , 

If component 1 is selected, draw X\ ^ TN^/ p+oo )(0,\) : X 2 \Xi = x\ ~ N(pXx, v 2 ) 1 
and accept simulated pair (xx, x 2 ) if x 2 > a 2 . Otherwise, draw X\ ~ TN[ ai}a2 / p } (9, v 2 ), 
and accept with probability 

and, upon acceptance, complete with 

X 2 |Xi =#1 - T^ a2i+oo) (pa;i,^ 2 ). 

We show formally in Appendix B2 that the acceptance rate of this algorithm is 
bounded from below by 0.22, and we found numerically that the optimal lower 
bound seems to be 1/2, see Section [ 
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3.5. Computational cost. The above algorithms, except algorithm S + , involve 
a few evaluations of function $, which is expensive. But such evaluations can 
be bypassed in most cases. In algorithm S~ for instance, given the expression of 
acceptance probability (|3.7[) . one should accept the proposed value x\ if and only if 



$ (rxi + s) > ut(xi) 

where u is an uniform variate, and the exact expression of r, s, and t are easily 
deduced from (|3.7[) . If good, fast approximations of $ are available, such that 
§L{x) < &(x) < $(aj), it is enough to check that §(rxi+s) > ut(xi) (resp. 
$ (rxi + s) < ut(xi)) to accept (resp. reject) x\. It is only when ut(x±) is very 
close to $ (rxi + s) that an exact evaluation of $ (rx\ + s) is required. 

Such fast, good approximations of $ may be deduced from the tables of our 
univariate algorithm, see Section 12.31 Specifically, and using the same notations 
as Section |2"UI let z = rx\ + s, and assume that z € [a m - m , a max ] , then one may 
set &(z) = A{jy z j h \ + 1), where A(i) > $(xj+i) denotes the total area of all the 
regions up to region i, and may be computed beforehand and hard-coded in the 
program like the other constants y^, di and so on. One may define similarly §i(z) 
using &(z) = 1 — <&(— z) . Wh en z £ [a max ,a m in], one may use the expansion of 
Abramowitz and Stegunl ( 1965 . p. 932) to derive the following upper and lower 



approximations 

z i z z z^ J z y z z z^ z° 

for z < 0. (For z > 0, similar formulae are obtained using $(z) = 1 — $(—«)). 
Note also that one does not necessarily have to compute all the terms of these 
expressions. For instance, if $>i(axi + 0) < uv(xi), where &i(z) = —<p(z)/z, then 
the proposed value can be rejected without computing the remaining terms. This 
principle can be used for each term of the expansion, but, on the other hand, it is 
preferable not to expand further the expressions above, since they work well already 
for reasonable values of a m ; n and a max , and since that would define diverging series. 

Provided the above strategy is implemented, the algorithms proposed in the 
section are reasonably fast, since they only involve a few basic operations, and 
their acceptance rate is greater than 1/2 for all parameters. Algorithms M + and 
M~ are slightly more expensive, as they require sampling a mixture index, but note 
that the same strategy can be implemented in order to avoid with good probability 
the evaluation of functions $ appearing in the expression of the mixture weights. 



3.6. Numerical illustration. We simulated 10 5 parameters (ai, a%, p), where p ~ 
U[— 1,1], ai, ai ~ N(0,s 2 ), conditional on ai > a^. For each parameter, we 
evaluated the acceptance rate of our algorithm by computing the average acceptance 
probability over the proposed draws generated by 1000 runs. Figure EP1 reports the 
histogram of the acceptance rates for s = 1; larger values of s give histograms that 
are even more concentrated around 1. 

In this simulation exercise, 90% of the acceptance rates are above 0.8 and 99% 
are above 0.65; none is lower than 1/2. 
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FIGURE 3.1. Histogram of acceptance rates corresponding to 10 5 
simulated vectors (a±,a2, p), where p ~ U[0, 1], a\, a,2 ~ iV(0, 1), 
conditional on a\ > a,2- 



4. Bl-DIMENSIONAL CASE: FINITE INTERVALS 

We now consider the simulation of X = (X 1: X 2 ) ~ N 2 {0, £), with 

1 P 
P 1 

conditional on X\ £ [ai,bi] and X2 € [02,62], where, without loss of generality, 
p > 0. As in the previous section, the main difficulty is to sample from the marginal 
density of X\ : 

(4.1) p(xi) oc tp{x{) [9(ax! + pi) - ®{axi + /3 )] J(oi < x x < h) 

where a = pjv > 0, fi\ = —a%/v, /3q = —bijv and /?i > /3o, since the distribution 
of X2 conditional on X\ = X\ is the univariate truncated Gaussian distribution 
TN [a2M] (pxi,v 2 ). 

We shall consider two situations, according to whether j3\ — /3q < A or not; we 
take A = 2, which seems to close to the optimal value in our simulations, in terms 
of minimum acceptance rate. 

When /?i — f3o is small (case T), one may Taylor expand the second factor of 
(|4.ip . which is denoted k from now on, into: 

k{xi) = §{axi + pi) - $(axi + p ) w {px - p )cp [a Xl + ^° | ^ 

hence k behaves like a Gaussian density, see the right panel of Figure POl Therefore, 
p(xi) is also well approximated by a Gaussian, which is the basis of algorithm T 
detailed in 4.2. 

Otherwise, when P\ — p is large, k behaves like the curve plotted in the left 
panel of Figure 14.11 In this case, called M 3 below, we cut K into at most three 
pieces, and derive a mixture proposal, using ideas similar to the previous section. 

4.1. Case A/ 3 . Let j t = —Pi/a, for i = 0, 1, and v = (Pi — fio)/2; note 71 < 70 
since Pq < pi. We may divide the curve of k into three parts, so as to re-use 
the same ideas as in Section [3l that is, deriving a mixture of Gaussian proposals. 
Specifically, 
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FIGURE 4.1. Function n{xi) — §{ax\ + ft) — <§>{ax\ + fti) , for 
(a, ft,, ft) = (2,-5,5) (left), and (a, A), ft) = (2,-0.5,0.5) (right) 



• for x\ < 71, one shows easily that 

{2$(A) - 1} 9{axi + ft) < k(xi) < $(ox x + ft). 

We know already that target density ip(x\)^(axi + px)I{xi < 71) can be 
simulated efficiently using a particular Gaussian proposal density, see Com- 
ponent 2 of Algorithm M + in Section f3T4l The above inequality indicates 
that, using the same proposal, one should obtain an acceptance rate which 
is at least 2$ (A) — 1 times the minimum acceptance rate of M + , that is, 
$(A) - 1/2 w 0.477. 

• for xi £ [71,70], K is roughly flat, and 

$(A) - 1/2 < k{xi) < 2$(u) - 1 < 1. 

This suggests using proposal distribution X\ ~ TiVr 7lVai)7oA ^ 1 i(0, 1), and 
accept realisation £1 with probability k{x\)/ {2$(-l>) — 1}. The acceptance 
probability is then bounded from below by 3>(A) — 1/2 0.477. 

• for x\ > 70, one has: 

{2$(A) - 1} $(-asi - ft) < k(3!i) < $(-axi - ft). 

Again, one may use the same proposal as for Component 2 of algorithm 
M~ , see !3.3| which should lead to an acceptance rate which is larger than 
$(A) - 1/2 w 0.477. 

The principle of Algorithm M 3 is therefore to draw from a mixture of at most three 
components, the relative weights of which are given below, and given the chosen 
component, to use one of the three strategies described above. 

Denote Q, £ C) £ r , the unnormalised weights of the left, centre, and right com- 
ponents, respectively. In case b\ V 71 > 0, one has 

f '1; 2V " 2 ) {* - * {^P) } m > «,) 







v d(a.a\ + ft) 
V5n 



exp 
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with mi = p{a,2 + \v) and function d was defined in Section I3.4( Otherwise, one 
obtains the same expression, but with A set to 0, i.e., 

^ exp (_|){ 4 (2^i)_ 4 (^)} /(7iaoi) 

with mi — pa2- In all cases, 

Cc = {2$(w) - 1} {$( 7o A bi) - $(71 V ai)} I{h > 71; a x < 70), 

One may show that the acceptance rate of this algorithm is bounded from below 
by 1/2 — <£>(— A) pa 0.477; simulations suggest this bound is optimal. We omit the 
exact calculations, as they are similar to those of previous algorithms. We managed 
to obtain this result under the following assumptions: (i) p > 0; (ii) b2 > and 
(iii) either 02 > a\ or b\ < 0. It is always possible to enforce such conditions, by 
either swapping X\ and X2, or changing their signs, or both. We note also that, in 
most of our simulated exercises, at least one component of this mixture is empty, 
and often two of them are, which makes it possible to skip calculating the weights 
and simulating the mixture index. 

4.2. Case T. As explained above, when f3\ — /3q < A, a good Gaussian approxim- 
ation of p(xi) is 



q(xi) oc ip(xi)ip \^axi 

that is, a N(m, s 2 ) density with 

/ 2n ( a(A)+/3i 
(to, s ) — ' 



2(1 + a 2 ) ' 1 + c 

In our experiments, this approximation appears to be accurate for all values of 
a, Po, P\ (in the sense that the rejection rate of the algorithm we now describe is 
always larger than 0.47 in our simulations, see next subsection). On the other hand, 
it seems difficult to apply approximations si milar to those we used before. Instead, 
we note that p{x\) is a log-concave density (jPrekopa . 1973). This suggests using 



either exponential or piecewise exponential pro posals, and working ou t a simplified 



version of ARS (Adaptive Rejection Sampling, iGilks and Wildl . 119921 ). 
Specifically, let £ denote the marginal log-density of X\ : 

£(xi) = k>g<p(xi) +logn(xx) 
the derivative of which is easy to compute: 

, _ _ Lp(ax 1 + ffi) - ipjax! + (3 ) 

^ {Xl> Xl+a $(ax 1 +p 1 )-<I>(ax 1 +M' 
which leads to the inequality 

S(xi)<S(v)+?(v)(xi-v) 

for any x\ , v £ [ai , 61] . Up to a constant, the right hand side is the log-density of the 
truncated Exponential distribution Exp[ Qi i, 1 i(A) defined in (|2.4p . with A = £'(u); 
note that A may be negative. Thus, one may sample x\ ~ Expj ai hi i {£'(w)}, and 
accept with probability 
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1.0 



Figure 4.2. Histogram of acceptance rates for 10 5 simulated para- 
meters (p, ai, b\, a2, 62), with p ~ U[—l, 1], 01, 02 ~ AT(0, 2 2 ), and 
&i = + 2e^ with &i ~ Exp(l), i = 1,2. 



Obviously, the difficulty is to choose v. Since the target density is well approximated 
by a TN\ ai fai (to, s 2 ), and assuming that a\>ra (resp. b\ < m), it seems reasonable 
to set v — (to + s) V ai (resp. v = (m — s) A 61). These values would be optimal 
if the target density would be equal to its approximation TWu^&ji (m, s 2 ). In case 
a\ < m and b\ > to, i.e., the mean to is within [ai,6i], we use instead a mixture 
proposal, based on two well chosen points v, w € [oj, 61], say v < w. Thus, 

£(zi) < + - w)} a UH + e'H(^i - w)} 

and one may sample from the density defined as the exponential of the right hand 
side above (which is a piecewise exponential distribution), and accept with prob- 
ability 

exp - {£(«) + e'^^i - «)} A {^H + e(w)( Xl - w)}] . 

Note that, in the original ARS algorithm o f lGilks and Wildl (|l992h . the proposal 



is progressively refined by adding a new component each time the proposed value 
is rejected. We found however that the acceptance rate of Algorithm T described 
above is generally above 1/2, so using fixed proposals with at most two components 
seems reasonable. Obviously, the good properties of the above algorithm lie in the 
good choice of points v, w, which was made possible by the knowledge of a good 
approximation of the target density. 

We were not able to prove formally that the acceptance rate of Algorithm T is 
bounded from below, as in previous cases, so we performed intensive simulations 
for assessing its properties; the acceptance rate of Algorithm T seems to converges 
to one for limiting values, say a± — > —00 but with bi — a\ and p kept fixed; and to 
be bounded from below by 1/2. 

4.3. Numerical illustration. We simulated 10 5 parameters (a%, b\, 0,2, 62, p), where 
p ~ U[-l, 1], a 1; a 2 ~ iV(0, 2 2 ), and b t = a* + 2a with e { - Exp(l), i = 1,2. For 
each parameter, we evaluated the acceptance rate of our algorithm by computing 
the average acceptance probability over the proposed draws generated by 1000 runs. 
Figure l4~2l reports the histogram of the acceptance rates. About 90% of these val- 
ues are above 0.71, about 99% are above 0.55, and all values are above 0.47, as 
expected. 
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5. Generalisation to 3 or more dimensions 

We discuss briefly the problem of simulating d-dimensional truncated Gaussian 
distributions, for d > 3 and for semi-finite truncation intervals; i.e., X ~ Nd(0d, X), 
subject to Xi > Oj, i = 1, . . . , d. As we have done for d = 2, we assume without 
loss of generality that a\ > . . . > a^,. It does not seem possible to generalise to 
dimension d > 3 algorithms based on mixture proposals, i.e., M + and M~ , as the 
corresponding mixture weights would involve intractable integrals. But algorithms 
S + and S~ can be generalised to larger dimensions, as explained below, which 
makes it possible to sample X under certain conditions on X and a = (oi, . . . , a<j). 

5.1. Extension of Algorithm S + . An obvious way of generalising Algorithm S + 
to 3 dimensions is to do the following. Let Q — (qij) — X -1 and X12 denote the 
sub-matrix obtained by removing the last row and the last column from X, then: 

(1) Sample {X\, X2) ~ ^(Cb, X12) conditional on Xi > a\ and X2 > a 2 (using 
one of the algorithms discussed in Section [3]) 

(2) Sample from the unconstrained conditional distribution of A3: 



X 3 \{X 1 =x 1 ,X 2 =x 2 }~ N 



913^1 + 923^2 1 
933 933 



(3) If X3 > 03, accept the simulated vector (xx, £2,^3); otherwise reject and go 
to Step 1. 

One easily shows that the acceptance rate corresponding to Step 3 is larger than 1/2 
under the following set of conditions: 913 < 0, 923 < and (71301 +92302 + 933^3 < 0. 

One may iterate the above principle so as to extend Algorithm S + to any di- 
mension d; i.e., for d = 4, add Step 4 where X4 is simulated from the appropriate 
conditional distribution and accept if X4 > 04; provided appropriate conditions 
similar to those above, are iteratively verified, one obtains an overall acceptance 
rate that is at least 2~( d ~ 1 \ since each rejection step (including Step 1 above for 
the first two variates X\ and X%) induces an acceptance rate that is at least 1/2. 

We note that these iterative conditions imply in particular that any pair of 
components of A is positively correlated, except for {X\, X%) which may have any 
type of correlation. There are several practical settings whe re this assumption is 



met, such as in Gaussian Markov random fields models (e.g. iRue and Heidi . 120051) 
where one would impose positive correlation between neighbour nodes. Thus, in 
such or other particular settings, and since the acceptance rate is expected to be 
higher than its lower bound for most parameters, as observed in two dimensions, 
the above algorithm may remain practical for dimensions as large as 4 or 5. For 
instance, in a Markov chain Monte Carlo context involving truncated Gaussian 
vectors or large dimension, one may try to form a larger and larger block, by 
including one variable at a time, checking the recursive assumptions above, and 
stop when either they are no longer met or the acceptance rate is too small. 

5.2. Extension of Algorithm S~ . Again, assuming d — 3, one notes that the 
marginal distribution of (Xi, X2) is: 



p(xi,x 2 ) cx exp --^ QijXiXj > $ 

{ i,j=l ) \ %S 



i3%i + 93303 



I (xi > ai;x 2 > a 2 ) 
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which suggests the following bivariate truncated Gaussian density as a proposal 
density: 



P 



j 1 ,A \ l^i=l 9i3^i T y33«3 I j 

(xi,x 2 ) oc exp <i -- QijXtXj Y q ) I ^ Xl - ai ' X2 - ° 2 ) 




based on inequality (|3.4[) . For given x\ and X2, the acceptance probability is there- 
fore: 

' "^ 2 ,i QiSXi + g 3 3Q3 \ , / Ya=1 laXi + 933«3 \ 

1/2 I / c I 1/2 I 

9 33 / V ?33 / 

where we recall that if)(x) — $(— x)/(p(x). Using the same type of calculations as in 
Section[331 one may show that the expectation of the acceptance probability above 
is larger than or equal to 1/2 provided q i3 > 0, (723 > 0, and 51301 +<723«2+933a3 > 0. 
Again, this means that the overall acceptance rate is larger than or equal to 1/4. 

As in the previous section, one may iterate the construction above, so as to 
obtain a simulation algorithm for any dimension d, the acceptance rate of which 
is bounded from below by 2~( d ~ 1 \ This requires checking recursively conditions 
similar to those above. The same remarks in the previous subsection relative to the 
applicability of this algorithm may be repeated here. 

6. Conclusion 

We focused in this paper on the simulation of independent truncated Gaus- 
sian variables, but similar ideas can be used in other settings, such as importance 
sampling or MCMC. For instance, in case T, see Section |4~2| one may use the de- 
rived Gaussian approximation as an importance distribution, rather than a basis 
of an ARS algorithm. The same remark applies to most of our algorithms. 

As briefly mentioned in the previous section, if one needs to simulate a vector 
from a high-dimensional truncated Gaussian distribution using MCMC, one may 
ask how to choose blocks of two or more variables, which will be updated using the 
algorithms proposed in this paper, in a way that ensures good MCMC convergence 
properties. A good strategy would be to form a first block of two variables with 
strong (conditional) correlation, then to see if additional variables may be included 
in that block, using the conditions given in the previous section, and repeat this 
process until all variables are included in a block of at least two variables. But more 
research is required to find the best trade-off in terms of convenience and efficiency. 

Open source programs implementing the proposed algorithms are available at 
the author's personal page on the website of his institution, www. crest .fr. 
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Appendix A: Outline of the univariate algorithm for a semi-finite 

TRUNCATION INTERVAL 

Note Devroye(a) refers to Devroye's algorithm, Direct (a) refers to the rejec- 
tion algorithm based on the non truncated Gaussian distribution, see Section ^. ll for 
details. In both cases the input a is the truncation point. Pre-computed constants 
consist of five floating-point tables: (xi), {yi), {y.), {di) and {Si); one integer table: 
(j'fc)i plus two design parameters a min , and a max . 



Require: a {truncation point} 
Ensure: x {simulated value} 
if a < a m i n then 

return Direct (a) 
else if a > a max then 
return Devroye(a) 
end if 

*a <~ 3\a/h\ 

loop 

Sample integer i uniformly between i a and N 
if i = N then {rightmost region} 

return Devroye(a^iv) 
else if i < i a + 1 then {two leftmost regions} 
Sample u - U[Q, 1] 
x — Xi + d{ * u 
if x > a then 

Sample v ~ U[0, 1] 

y <- yi * v 

if y <y. then 

return x 
else if y < tp(x) then 

return x 
end if 
end if 

else {all the other regions} 
Sample u - U[0, 1] 

y <- u * yi 

if y < y ■ then {occurs with high probability} 

return Xi + u * Si 
else 

Sample v - U[0, 1] 
x <— Xi + di * v 
if y < (p(x) then 

return x 
end if 
end if 
end if 
end loop 
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Appendix B: Lower bounds for Acceptance rates 

Bl. algorithm S~ 

Let A(a\,a 2 , p) the acceptance rate (|3.8I) . which we rewrite as: 

A(a 1 ,a 2 ,p) = £tat [Qi=o) (^ p 2) bP {Z)\ /c (-a) 

where Z = —{pX\ — a 2 )/v ~ TN\ a oo -dfi, p 2 ), a = (a 2 — pa\)/v, and ft = a 2 v. Note 
that a > 0, ft < a, and tf> is a decreasing function. Thus, the quantity above is a 
decreasing function of ft. (To see this, one can rewrite the distribution of Z as 

Z = ft + p®- 1 (u + (i- ^) 



where U is an uniform variate, and check that, conditional on U — u, Z is a 
decreasing function of ft.). Thus, the above quantity is larger than or equal to the 
same quantity, but with ft = a: 

A(a 1 ,a 2 ,p) > E T N [a}00) ( a , p 2) [i/j(Z)] /c(-a) = E TN[0aa)(0A) [ip (pZ' + a)] /c(-a) 

which is a decreasing function of p, hence 

(6-1) A(a 1 ,a 2 ,p) > £W [0i0o) (o,i) [ip (Z' + a)} /c(-a) , 

and since c(— a) = (yjit /2) A (1/a), one can show that the bound is minimised for 
a = which leads to: 



A(ai,a 2 ,p) > \ -E TN 



)(0,1) 



0.416. 



This lower bound is not sharp, because not all combinations of (a, ft, p) are valid, 
even in the constraints a > 0, ft < a are taken into account; for instance, ft = a 
implies that <zi = pa 2 < p 2 a\, which is impossible if p ^ 1. Our simulations suggests 
that the optimal lower bound is 1/2, see Section [3T6l 

B2. ALGORITHM M + 

One easily shows that x( x )/x( x ') > 1/2 for all x,x' £ [0,Xd], with Xd ~ 3.117. 
Thus, if (a 2 — pa{)/v < Xd, the acceptance rate is larger than or equal to 1/2 by 
construction. Now assume that (a 2 — pa\)/v > Xd', note that x( x ) 1S an increasing 
function for x > Xd- Since a\ > a 2 and p < 1, one has 

9 = p(a 2 + Xv) < a 2 + v{\p — Xd) < a\. 

The acceptance rate equals 



(6.2) 



3 



= E- 



X(Z) 



X (Zmax) 



(a 2 - pai)/v, and r) = (a 2 - p0)/u > z max ; note d(z max ) = x( z max) 



> 0.751, but we assumed that z mnx > Xd ~ 3.117. The TN\Q^ a ^_ pai y v ^ 



where z max 

provided z max ^> u. 101, out we assumed mat z max ,> Xd ~ o.n i . uiei JV [o,(a 2 -pa 
distribution should concentrate its mass at the right edge of interval [0,z max ], 
and x(Z)/x(z max ) should take values close to one. Specifically, one has that 
z$(-z)/ip(z) E (0.84, 1) for z > 2, thus 



X^max) 



> 



0.84e A( ^ Zmax) > 0.5 
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for z E [zmax — 0.76, z max ]- Therefore (|6.2II is larger than 0.5 times the probability 
that Z > Zmax — 0.76, for Z ~ TN\q Zx i (77, p 2 ), which is larger than or equal to 0.44, 
the probability of the same event with respect to Z ~ TNu3 tZmxx ](z max , 1), for z max . 
This gives a lower bound for (|6.2p of 0.22. We obtained sharper bounds with more 
tedious calculations (omitted here), but more importantly, our simulation studies 
indicates that the optimal lower bound is likely to be larger than or equal to 1/2. 



